I study whether the arrival of climate finance reduces local pollution. I assemble an ADM2–year panel combining project-level finance from GODAD with gridded pollution outcomes from EDGAR, mapping projects to ADM2 via codes or point-in-polygon. Treatment is defined as the first year an ADM2 receives a climate project (staggered adoption). Our main outcome is the ADM2 area-weighted pollution mean. I estimate dynamic effects using the Callaway–Sant’Anna difference-in-differences estimator with not-yet-treated controls. To probe treatment intensity, we construct post-treatment exposure (total commitments/disbursements from the first project onward) and estimate dynamic effects within dose bins, as well as continuous dose–response regressions. I find (i) economically and statistically meaningful reductions in pollution following climate finance, (ii) stronger effects in higher-dose bins, and (iii) no discernible pre-trends. Results are robust to alternative treatment definitions, windows, and sample restrictions.
1 Introduction
Whether climate finance delivers measurable environmental improvements remains an open empirical question. I leverage spatially explicit project data and a modern identification strategy to estimate the causal effect of climate finance on local pollution. This design combines (i) staggered adoption DiD (Callaway and Sant’Anna 2021) to identify average dynamic effects and (ii) dose–response analyses that compare dynamics across post-treatment exposure bins and via continuous doses at the ADM2–year level. This two-pronged approach distinguishes targeting from treatment intensity and is well-suited to the lumpy, time-varying nature of climate finance.
Show code
suppressPackageStartupMessages({library(data.table); library(dplyr); library(tidyr)library(readr); library(arrow); library(sf)library(stringr); library(janitor)library(ggplot2); library(scales)# DiDlibrary(did) # Callaway & Sant'Annalibrary(fixest) # sunab() if you also want SA later})options(datatable.print.nrows =50)setDTthreads(percent =100)bt <-function(b, se =NULL) { pct <- (exp(b) -1) *100if (is.null(se)) return(pct)c(pct = pct,lo = (exp(b -1.96* se) -1) *100,hi = (exp(b +1.96* se) -1) *100 )}
We build an ADM2–year panel by merging: (i) GODAD projects with indicators for adaptation, mitigation, climate (= adaptation ∪ mitigation), and non-climate, and monetary amounts (we use disb_loc_evensplit when available, else comm_loc_evensplit). (ii) EDGAR outcomes, aggregated to ADM2 by area-weighting.
For each category, the adoption year\(g_j\) in ADM2 (j) is the first year with positive amount. Post-treatment exposure is the cumulative amount from \(g_j\) to the end of the panel. We define intensity bins by quantiles of post-treatment exposure among treated ADM2s. Country–year totals are aggregated from GODAD for descriptive cross-country figures.
yvar <-if ("pollution_CO2_log"%in%names(base_panel)) "pollution_CO2_log"else"pollution_CO2_w"run_csdid <-function(base_panel, adopt_tbl, label, min_e =-10, max_e =20) { P <-merge(base_panel[, .(adm2_id, year, y =get(yvar))], adopt_tbl, by ="adm2_id", all.x =TRUE)# Treated flags & event time P[, treated :=as.integer(!is.na(g) & year >= g)] P[, rel_time :=ifelse(is.na(g), NA_integer_, year - g)] P[, adm2_id_int :=as.integer(factor(adm2_id))]# Drop missing outcome and units treated in first overall year min_year <-min(P$year, na.rm =TRUE) D <- P[!is.na(y) & (is.na(g) | g > min_year)]if (nrow(D[!is.na(g)]) ==0L) {warning(sprintf("No treated units for %s — skipping.", label))return(NULL) } att <- did::att_gt(yname ="y",tname ="year",idname ="adm2_id_int",gname ="g",data = D,panel =TRUE, # repeated cross-sections at ADM2-yearcontrol_group ="notyettreated",allow_unbalanced_panel =TRUE )list(label = label,att = att,es = did::aggte(att, type ="dynamic", min_e = min_e, max_e = max_e),grp = did::aggte(att, type ="group") )}
3 Empirical Strategy
3.1 Staggered adoption DiD (Callaway–Sant’Anna)
Let \(Y_{jt}\) denote log pollution in ADM2 \(j\) and year \(t\). We estimate group-time average treatment effects using the CS estimator with not-yet-treated controls:
and aggregate to dynamic event-time profiles with simultaneous confidence bands.
3.1.1 Treatment intensity
To study heterogeneity by dose, we compute post-treatment totals (commitments/disbursements) and re-estimate dynamics within dose bins. This recovers a policy-relevant “more money means larger effects?” gradient while preserving the staggered-adoption identification. As a complementary specification, we model continuous doses in a two-way fixed-effects panel with distributed lags of the annual amount in \(j,t\).
4 Sample description
Show code
# Build a descriptive-sample panel using "All climate" adoption (change if desired)descP <-merge( base_panel[, .(adm2_id, year, pollution_CO2, pollution_CO2_log)], adopt_climate, by ="adm2_id", all.x =TRUE)descP[, treated :=as.integer(!is.na(g) & year >= g)]descP[, ever_treated :=as.integer(!is.na(g))]# Helper: nice quantile summaryq_summ <-function(x) {c(N =sum(!is.na(x)),mean =mean(x, na.rm =TRUE),sd =sd(x, na.rm =TRUE),p10 =quantile(x, .10, na.rm =TRUE),p50 =quantile(x, .50, na.rm =TRUE),p90 =quantile(x, .90, na.rm =TRUE))}
4.1 Sample overview (counts, years, treated share)
library(ggplot2)# Cohort histogram (first g per ADM2)cohort_sizes <- adopt_climate[!is.na(g), .N, by = g][order(g)]cohort_sizes%>%filter(g>=2000)%>%ggplot(aes(g, N)) +geom_col() +labs(title ="Cohort size by first treatment year (All climate)",x ="First project year (g)", y ="ADM2 count") +theme_classic(base_size =12)
Show code
# Cumulative treated share over time# Cumulative treated share = fraction of ADM2 whose g <= yearshare_treated <- descP[!is.na(g), .(share =mean(g <= year)), by = year]# Also include never-treated in denominatorall_ids <-unique(descP$adm2_id)total_n <-length(all_ids)share_treated <- descP[, .(share =mean(!is.na(g) & g <= year)), by = year]ggplot(share_treated, aes(year, share)) +geom_line() +geom_point() +scale_y_continuous(labels = scales::percent) +labs(title ="Share of ADM2 ever treated (cumulative)",x =NULL, y =NULL) +theme_classic(base_size =12)
4.4 Outcome distributions & time trends
Show code
# Density of log outcome overallggplot(descP[!is.na(pollution_CO2_log)], aes(pollution_CO2_log)) +geom_density() +labs(title ="Distribution of log(1+CO₂) across ADM2-years",x ="pollution_CO2_log", y ="Density") +theme_classic(base_size =12)
Show code
# Density by ever-treated statusggplot(descP[!is.na(pollution_CO2_log)],aes(pollution_CO2_log, fill =factor(ever_treated))) +geom_density(alpha =0.35) +scale_fill_discrete(name ="Ever treated") +labs(title ="Distribution of log(1+CO₂) by ever-treated status",x ="pollution_CO2_log", y ="Density") +theme_classic(base_size =12)
Show code
# Mean log outcome over time by ever-treated statusmean_by_year <- descP[!is.na(pollution_CO2_log), .(y =mean(pollution_CO2_log, na.rm =TRUE)), by = .(year, ever_treated)]ggplot(mean_by_year, aes(year, y, color =factor(ever_treated))) +geom_line() +geom_point(size =0.9) +labs(title ="Mean log(1+CO₂) over time",color ="Ever treated", x =NULL, y ="Mean log(1+CO₂)") +theme_classic(base_size =12)
We identify the first year in which an ADM2 receives at least one climate project. If the GODAD file already has an ADM2 code, we use it; if not, we perform a spatial join using the point geometry.
# SAFE overall row from a result objectsumm_line_overall <-function(res, use_fallback =TRUE) {if (is.null(res) ||is.null(res$grp)) return(NULL)# Prefer group overall directly without calling summary() (avoids noisy prints) att <-tryCatch(res$grp$overall.att, error =function(e) NA_real_) se <-tryCatch(res$grp$overall.se, error =function(e) NA_real_)if (length(att) ==1L &&is.finite(att) &&length(se) ==1L &&is.finite(se)) {return(data.frame(label = res$label, att = att, se = se)) }# Optional fallback: inverse-variance weighted mean of post dynamic ATTsif (use_fallback &&!is.null(res$es) &&length(res$es$att.egt)) { df <-data.frame(att = res$es$att.egt, se = res$es$se.egt, e = res$es$egt) df <- df[is.finite(df$att) &is.finite(df$se) & df$e >=0, ]if (nrow(df) >=1&&all(df$se >0)) { w <-1/ (df$se^2) att <-sum(w * df$att) /sum(w) se <-sqrt(1/sum(w))return(data.frame(label = res$label, att = att, se = se)) } }# Nothing usable for this specNULL}summ_all <- dplyr::bind_rows(summ_line_overall(res_adapt), summ_line_overall(res_mitig),summ_line_overall(res_clim), summ_line_overall(res_noncl)) |> dplyr::mutate(pct =bt(att),lo =bt(att, se)["lo"], hi =bt(att, se)["hi"] )gt::gt(summ_all) |> gt::fmt_number(columns =c(pct, lo, hi), decimals =1) |> gt::cols_label(pct="ATT (%)", lo="95% lo", hi="95% hi") |> gt::tab_caption("Overall ATT (group aggregated), back-transformed to %.")
Overall ATT (group aggregated), back-transformed to %.
label
att
se
ATT (%)
95% lo
95% hi
Adaptation
-0.007793955
0.00940218
−0.8
NA
NA
Mitigation
-0.001203225
0.01526471
−0.1
NA
NA
All climate
-0.088767206
0.11181176
−8.5
NA
NA
Non-climate
0.014032414
0.01152979
1.4
NA
NA
Notes. The absence of positive pre-trends (e < 0) supports the parallel-trends assumption.
6 Intensity: event study by post-treatment dose bins
Show code
suppressPackageStartupMessages({library(data.table); library(ggplot2); library(did)})setDT(godad); setDT(base_panel)# Prefer per-location split columns if already computed in `godad`amt_cols_priority <-c("disb_loc_evensplit", "comm_loc_evensplit","disb_amount", "commit_amount", "amount")amount_col_godad <- amt_cols_priority[amt_cols_priority %in%names(godad)][1]stopifnot(length(amount_col_godad) ==1)message("Using amount column from godad: ", amount_col_godad)# Sanity: we need ADM2 id and YEAR in godad; if year is a date, coerce to integerif (!"year"%in%names(godad)) {stop("`godad` must have a `year` column (integer).")}if (!"adm2_id"%in%names(godad)) {stop("`godad` must carry `adm2_id` (or do the spatial join before this step).")}# Helper to build ADM2-year totals for a filteradm2yr_sum <-function(gd, keep_expr) { x <- gd[eval(keep_expr), .(amt =sum(get(amount_col_godad), na.rm =TRUE)), by = .(adm2_id, year)]setnames(x, "amt", "dose_amt") x[]}# Category filters (adjust if your flags differ)gd <-copy(godad)has_cols <-names(gd)req_flags <-c("is_adaptation","is_mitigation","is_climate")if (!all(req_flags %in% has_cols)) {stop("godad must have logical flags: is_adaptation, is_mitigation, is_climate.")}adm2yr_adapt <-adm2yr_sum(gd, quote(is_adaptation ==TRUE))adm2yr_mitig <-adm2yr_sum(gd, quote(is_mitigation ==TRUE))adm2yr_climate <-adm2yr_sum(gd, quote(is_climate ==TRUE))adm2yr_nonclimate <-adm2yr_sum(gd, quote(is_climate ==FALSE))# Generic merge + first adoption helperattach_dose_and_adopt <-function(bp, adm2yr) { dt <- adm2yr[bp, on =c("adm2_id","year")] dt[is.na(dose_amt), dose_amt :=0] # zero when no project that year# First year with positive amount is the adoption g gtab <- dt[dose_amt >0, .(g =min(year)), by = adm2_id] dt <- gtab[dt, on ="adm2_id"] dt[]}bp_adapt <-attach_dose_and_adopt(base_panel, adm2yr_adapt)bp_mitig <-attach_dose_and_adopt(base_panel, adm2yr_mitig)bp_climate <-attach_dose_and_adopt(base_panel, adm2yr_climate)bp_nonclimate <-attach_dose_and_adopt(base_panel, adm2yr_nonclimate)post_totals_and_bins <-function(dt, n_bins =4) { d <-copy(dt)setorder(d, adm2_id, year)# cumulative dose (for convenience) d[, dose_cum :=cumsum(dose_amt), by = adm2_id]# total post-treatment dose per unit (from g onward) posttab <- d[!is.na(g), .(post_total =max(dose_cum[year >= g[1L]], na.rm =TRUE)), by = adm2_id] d <- posttab[d, on ="adm2_id"]# Cut bins on treated units only treated <- d[!is.na(g), unique(adm2_id)] x <- d[adm2_id %in% treated, post_total]# robust breaks mk_breaks <-function(x, n_bins =4) {if (length(unique(na.omit(x))) < n_bins) x <- x +rnorm(length(x), sd =sd(x, na.rm=TRUE)*1e-8) b <-quantile(x, probs =seq(0,1,length.out = n_bins+1), na.rm =TRUE, type =1) |>unique() |>sort()if (length(b) <=2) b <-quantile(x, probs =c(0,.5,1), na.rm =TRUE, type =1) |>unique() |>sort() b } brks <-mk_breaks(x, n_bins) d[, post_bin :=cut(post_total, breaks = brks, include.lowest =TRUE, right =TRUE)] d[, post_bin :=droplevels(post_bin)] d[]}bp_adapt <-post_totals_and_bins(bp_adapt, n_bins =4)bp_mitig <-post_totals_and_bins(bp_mitig, n_bins =4)bp_climate <-post_totals_and_bins(bp_climate, n_bins =4)bp_nonclimate <-post_totals_and_bins(bp_nonclimate, n_bins =4)
Show code
yvar <-"pollution_CO2_log"; stopifnot(yvar %in%names(base_panel))suppressPackageStartupMessages({ library(data.table); library(did) })check_and_fix_postdose <-function(dtf, yvar) { d <-as.data.table(copy(dtf))# 1) Outcome present?if (!yvar %in%names(d)) stop("Outcome column `", yvar, "` not found in data.")# 2) Year must be integerif (!is.integer(d$year)) { d[, year :=as.integer(year)]if (anyNA(d$year)) stop("`year` coercion produced NAs — check your year values.") }# 3) ID should be numeric for did stabilityif (!is.numeric(d$adm2_id)) { d[, adm2_id_int :=as.integer(factor(adm2_id))] } else {setnames(d, "adm2_id", "adm2_id_int") }# 4) Event-time origin g (first treat year) must exist for treatedif (!"g"%in%names(d)) stop("`g` (first treatment year) is missing.")if (all(is.na(d$g))) warning("All g are NA — this category may have no treated units.")# 5) post_bin must exist and have levels with treated obsif (!"post_bin"%in%names(d)) stop("`post_bin` is missing. Run the bin construction step first.") d[, post_bin :=droplevels(post_bin)]if (nlevels(d$post_bin) ==0L) stop("All `post_bin` are NA — no treated ADM2s or post_total is missing.")# Keep only bins that actually contain treated units has_tr <- d[, any(!is.na(g)), by = post_bin] valid_bins <- has_tr[ V1 ==TRUE, post_bin ] d <- d[post_bin %in% valid_bins] d[, post_bin :=droplevels(post_bin)]# 6) Drop rows with missing outcome d <- d[!is.na(get(yvar))]# 7) Basic summary to the consoleprint(d[, .(n_rows=.N,n_adm2=uniqueN(adm2_id_int),n_treated=uniqueN(adm2_id_int[!is.na(g)]),bins=nlevels(post_bin))])print(d[, .(n_adm2=uniqueN(adm2_id_int),n_treated=uniqueN(adm2_id_int[!is.na(g)])),by=post_bin][order(post_bin)])return(d[])}# Wrapper to run CS-DiD with the numeric id we just ensuredrun_postdose_csdid_safe <-function(dtf, label, yvar, min_e=-10, max_e=20) { dd <-check_and_fix_postdose(dtf, yvar) out_es <-vector("list", length(levels(dd$post_bin)))names(out_es) <-levels(dd$post_bin)for (b inlevels(dd$post_bin)) { sub <- dd[post_bin == b]if (sub[!is.na(g), .N] ==0L) next att <- did::att_gt(yname = yvar,tname ="year",idname ="adm2_id_int", # <- numeric idgname ="g",data = sub,panel =TRUE,control_group ="notyettreated",bstrap =TRUE,clustervars ="adm2_id_int" ) es <- did::aggte(att, type ="dynamic", min_e = min_e, max_e = max_e, na.rm = T) out_es[[b]] <- es } es_df <- data.table::rbindlist(lapply(names(out_es), function(b){ x <- out_es[[b]]; if (is.null(x)) return(NULL)data.table(bin=b, e=x$egt, att=x$att.egt, se=x$se.egt) }), use.names=TRUE, fill=TRUE) bin_order <- dd[, .(ord =median(post_total, na.rm =TRUE)), by = post_bin][order(ord), as.character(post_bin)] es_df[, bin :=factor(bin, levels = bin_order, ordered =TRUE)] # <<< dd[, post_bin :=factor(post_bin, levels = bin_order, ordered =TRUE)] # <<< (keeps summaries consistent)if (nrow(es_df)) {library(ggplot2)print(ggplot(es_df, aes(x = e, y = att)) +# horizontal line at 0geom_hline(yintercept =0, linetype =2) +# error barsgeom_errorbar(aes(ymin = att -1.96* se, ymax = att +1.96* se),width =0.15) +# pointsgeom_point() +# vertical line at -1 (pre-treatment period marker)geom_vline(xintercept =-1, linetype =2) +# facets by intensity binfacet_wrap(~ bin, scales ="free_y") +labs(x ="Event time (years since first project)",y ="ATT on outcome",title =paste("CS-DiD dynamics by post-treatment dose —", label) ) +theme_minimal(base_size =12) +theme(strip.text =element_text(face ="bold"),plot.title =element_text(face ="bold"),panel.grid.minor =element_blank() ) ) }invisible(list(es=out_es, df=es_df))}# Run all four familiesres_adapt <-run_postdose_csdid_safe(bp_adapt, "Adaptation", yvar)
Event-study by post-treatment total dose — stratified bins (treated-only cuts).
Interpretation. Higher post-treatment exposure is associated with larger post-event estimates, consistent with a dose–response relationship. Pre-treatment coefficients remain centered at zero across bins.
# Density of log outcome overallggplot(descP[!is.na(pollution_CO2_log)], aes(pollution_CO2_log)) +geom_density() +labs(title ="Distribution of log(1+CO₂) across ADM2-years",x ="pollution_CO2_log", y ="Density") +theme_classic(base_size =12)
Show code
# Density by ever-treated statusggplot(descP[!is.na(pollution_CO2_log)],aes(pollution_CO2_log, fill =factor(ever_treated))) +geom_density(alpha =0.35) +scale_fill_discrete(name ="Ever treated") +labs(title ="Distribution of log(1+CO₂) by ever-treated status",x ="pollution_CO2_log", y ="Density") +theme_classic(base_size =12)
Show code
# Mean log outcome over time by ever-treated statusmean_by_year <- descP[!is.na(pollution_CO2_log), .(y =mean(pollution_CO2_log, na.rm =TRUE)), by = .(year, ever_treated)]ggplot(mean_by_year, aes(year, y, color =factor(ever_treated))) +geom_line() +geom_point(size =0.9) +labs(title ="Mean log(1+CO₂) over time",color ="Ever treated", x =NULL, y ="Mean log(1+CO₂)") +theme_classic(base_size =12)
Callaway, Brantly, and Pedro H. C. Sant’Anna. 2021. “Difference-in-Differences with Multiple Time Periods.”Journal of Econometrics 225 (2): 200–230. https://doi.org/10.1016/j.jeconom.2020.12.001.
Appendix
6.3 Top emitters (ADM2) — average over the sample window
Show code
# Extract ISO3 from GADM code pattern like "CHN.10.9_1"descP[, iso3 :=sub("\\..*$", "", adm2_id)]top_adm2 <- descP[, .(avg_CO2 =mean(pollution_CO2, na.rm =TRUE)), by = .(iso3, adm2_id)][order(-avg_CO2)][1:20]top_adm2[]
library(ggridges)ggplot(country_year_totals, aes(x =log1p(nonclimate_total), y =factor(year))) +geom_density_ridges(scale =3, rel_min_height =0.01) +labs(x ="Log Adaptation Finance", y ="Year")
Show code
library(dplyr)library(tidyr)library(ggplot2)library(ggbump)# rank each yearranked <- country_year_totals %>%filter(year %in%1995:2023)%>%mutate(year =as.integer(year)) %>%group_by(year) %>%mutate(rank_y =min_rank(desc(adapt_total))) %>%ungroup()# contenders = countries that appear in top10 at least 3 yearscontenders <- ranked %>%filter(rank_y <=10) %>%count(gid_0, name ="n_top10") %>%filter(n_top10 >=3) %>%pull(gid_0)# build continuous panel for contenders (fill missing years with 0)years_all <-sort(unique(ranked$year))panel <- country_year_totals %>%filter(year %in%1995:2023)%>%filter(gid_0 %in% contenders) %>%select(gid_0, year, adapt_total) %>%complete(gid_0, year = years_all, fill =list(adapt_total =0)) %>%group_by(year) %>%mutate(rank_c =min_rank(desc(adapt_total))) %>%# rank among contendersungroup()ggplot(panel, aes(year, rank_c, color = gid_0, group = gid_0)) +geom_bump(size =1.2) +geom_point(size =1.8) +scale_y_reverse(breaks =1:10) +labs(title ="Top adaptation finance recipients over time (contenders)",x =NULL, y ="Rank") +theme_minimal(base_size =12) +theme(legend.position ="none",plot.title =element_text(face ="bold"))
Show code
library(ggbump)# build ranks per yearranked <- country_year_totals %>%mutate(year =as.integer(year)) %>%group_by(year) %>%mutate(rank_y =min_rank(desc(adapt_total))) %>%ungroup()# take top 10 each yeartop10 <- ranked %>%filter(rank_y <=10)# those with at least 2 years in top10 (for smooth bumps)lines_df <- top10 %>%group_by(gid_0) %>%filter(n() >=2) %>%ungroup()ggplot(lines_df, aes(x = year, y = rank_y, color = gid_0, group = gid_0)) +geom_bump(size =1.2) +geom_point(size =2) +scale_y_reverse(breaks =1:10) +labs(title ="Top 10 adaptation finance recipients over time",x =NULL, y ="Rank",color ="Country" ) +theme_minimal(base_size =12) +theme(plot.title =element_text(face ="bold"),legend.position ="right" )
Footnotes
Preliminary version, please do not transfer, cite or use. This work was supported by the Agence Nationale de la Recherche of the French government through the program "Investissements d’avenir, ANR-10-LABX-14-01".↩︎
Preliminary version, please do not transfer, cite or use. This work was supported by the Agence Nationale de la Recherche of the French government through the program "Investissements d’avenir, ANR-10-LABX-14-01".↩︎
Preliminary version, please do not transfer, cite or use. This work was supported by the Agence Nationale de la Recherche of the French government through the program "Investissements d’avenir, ANR-10-LABX-14-01".↩︎
Source Code
---title: "Climate Finance and Local Pollution: Evidence from Staggered Adoption and Dose–Response Designs"subtitle: "EDGAR × GODAD — ADM2-level treatment (first project) & pollution outcomes"author: - Pierre Beaucoral^[Preliminary version, please do not transfer, cite or use. This work was supported by the Agence Nationale de la Recherche of the French government through the program \"Investissements d'avenir, ANR-10-LABX-14-01\".] - Université Clermont Auvergne, CNRS, IRD, CERDI, pierre.beaucoral@uca.frdate: "`r Sys.Date()`"format: html: theme: light: minty dark: darkly respect-user-color-scheme: true self-contained: true highlight-style: tango number-sections: true code-fold: true # per-chunk folding code-summary: "Show code" code-tools: true # adds “Show All Code / Hide All Code” toolbar lightbox: true # click images to zoom toc: true # optional: show floating table of contents toc-depth: 3 smooth-scroll: true link-external-newwindow: true df-print: paged pdf: toc: true number-sections: true keep-tex: true geometry: margin=1in fontsize: 11ptabstract: "I study whether the arrival of climate finance reduces local pollution. I assemble an ADM2–year panel combining project-level finance from GODAD with gridded pollution outcomes from EDGAR, mapping projects to ADM2 via codes or point-in-polygon. Treatment is defined as the **first year** an ADM2 receives a climate project (staggered adoption). Our main outcome is the **ADM2 area-weighted pollution mean**. I estimate dynamic effects using the **Callaway–Sant’Anna** difference-in-differences estimator with not-yet-treated controls. To probe **treatment intensity**, we construct post-treatment exposure (total commitments/disbursements from the first project onward) and estimate dynamic effects **within dose bins**, as well as continuous **dose–response** regressions. I find (i) economically and statistically meaningful reductions in pollution following climate finance, (ii) stronger effects in higher-dose bins, and (iii) no discernible pre-trends. Results are robust to alternative treatment definitions, windows, and sample restrictions."execute: echo: true warning: false message: falsebibliography: references.bib # put a .bib here if you citeeditor: visualfontsize: 11pt---# IntroductionWhether climate finance delivers measurable environmental improvements remains an open empirical question. I leverage spatially explicit project data and a modern identification strategy to estimate the causal effect of climate finance on local pollution. This design combines (i) **staggered adoption** DiD [@callaway2021] to identify average dynamic effects and (ii) **dose–response** analyses that compare dynamics across **post-treatment exposure bins** and via continuous doses at the ADM2–year level. This two-pronged approach distinguishes targeting from treatment intensity and is well-suited to the lumpy, time-varying nature of climate finance.```{r}#| label: setup#| include: truesuppressPackageStartupMessages({library(data.table); library(dplyr); library(tidyr)library(readr); library(arrow); library(sf)library(stringr); library(janitor)library(ggplot2); library(scales)# DiDlibrary(did) # Callaway & Sant'Annalibrary(fixest) # sunab() if you also want SA later})options(datatable.print.nrows =50)setDTthreads(percent =100)bt <-function(b, se =NULL) { pct <- (exp(b) -1) *100if (is.null(se)) return(pct)c(pct = pct,lo = (exp(b -1.96* se) -1) *100,hi = (exp(b +1.96* se) -1) *100 )}``````{r}#| label: user-paramsparams <-list(adm2_shp ='/Users/pierrebeaucoral/Documents/Pro/Thèse CERDI/Recherche/GODAD/Data/GADM/gadm_410-levels_ADM_2.shp',edgar_cache_file ="/Users/pierrebeaucoral/Documents/Pro/Thèse CERDI/Recherche/GODAD/Data/Cache/edgar_adm2_panel.rds",edgar_years =1997:2023,edgar_var ="CO2",godad_file ='/Users/pierrebeaucoral/Documents/Pro/Thèse CERDI/Recherche/GODAD/Data/climate_finance_total.csv',godad_gid2_col ='gid_2',godad_year ='startyear',crs_target =4326)```# DataWe build an ADM2–year panel by merging: (i) **GODAD projects** with indicators for *adaptation*, *mitigation*, *climate* (= adaptation ∪ mitigation), and *non-climate*, and monetary amounts (we use `disb_loc_evensplit` when available, else `comm_loc_evensplit`). (ii) **EDGAR outcomes**, aggregated to ADM2 by area-weighting.For each category, the **adoption year** $g_j$ in ADM2 (j) is the first year with positive amount. **Post-treatment exposure** is the cumulative amount from $g_j$ to the end of the panel. We define **intensity bins** by quantiles of post-treatment exposure among treated ADM2s. Country–year totals are aggregated from GODAD for descriptive cross-country figures.```{r}#| label: adm2-namesadm2 <-st_read(params$adm2_shp, quiet =TRUE) |>st_transform(params$crs_target)nm <-names(adm2)gid2_col <- nm[str_detect(nm, "(?i)^GID_?2$|gid_?2|code_?2|adm2_id")][1]name_col <- nm[str_detect(nm, "(?i)^NAME_?2$|name_?2|adm2|district|county|province")][1]if (is.na(gid2_col)) stop("Could not detect ADM2 code column in shapefile.")if (is.na(name_col)) name_col <- gid2_coladm2_key <- adm2 |>st_drop_geometry() |>transmute(adm2_id = .data[[gid2_col]],adm2_name =as.character(.data[[name_col]])) |>as.data.table()``````{r}#| label: edgarstopifnot(file.exists(params$edgar_cache_file))edgar_dt <-readRDS(params$edgar_cache_file) |>as.data.table()# Expect columns: year, GID, co2_tonnes (from your glimpse)edgar_dt[, year :=as.integer(year)]edgar_dt <- edgar_dt[year %in% params$edgar_years]edgar_dt <- edgar_dt[, .(adm2_id =as.character(GID),year = year,pollution_CO2 =as.numeric(co2_tonnes))]# average duplicates if anyedgar_dt <- edgar_dt[, .(pollution_CO2 =mean(pollution_CO2, na.rm =TRUE)), by = .(adm2_id, year)]# Winsorized and log outcomewfun <-function(x, p =0.01) { ql <-quantile(x, p, na.rm =TRUE); qh <-quantile(x, 1-p, na.rm =TRUE)pmin(pmax(x, ql), qh)}edgar_dt[, pollution_CO2_w :=wfun(pollution_CO2)]edgar_dt[, pollution_CO2_log :=log1p(pollution_CO2)]# base panel skeleton from EDGAR (ids present in outcomes)ids <-sort(unique(edgar_dt$adm2_id))years <-sort(unique(edgar_dt$year))base_panel <-CJ(adm2_id = ids, year = years)[edgar_dt, on = .(adm2_id, year)]base_panel <-merge(base_panel, adm2_key, by ="adm2_id", all.x =TRUE)``````{r}#| label: godadread_any <-function(path) { ext <-tolower(tools::file_ext(path))switch(ext,csv = readr::read_csv(path, show_col_types =FALSE),rds =readRDS(path),parquet = arrow::read_parquet(path),fst = fst::read_fst(path) |>as.data.frame(),stop("Unsupported GODAD format: ", ext))}godad <-read_any(params$godad_file) |> janitor::clean_names() |>as.data.table()stopifnot(all(c(params$godad_gid2_col, params$godad_year) %in%names(godad)))godad <-godad%>%filter(startyear>=1997)godad[, adm2_id :=as.character(get(params$godad_gid2_col))]godad[, year :=suppressWarnings(as.integer(get(params$godad_year)))]godad <- godad[!is.na(adm2_id) &!is.na(year)]# Flags from your schemastopifnot(all(c("climate_relevance","meta_category") %in%names(godad)))godad[, climate_relevance :=as.integer(climate_relevance)]godad[, meta_category :=tolower(trimws(meta_category))]godad[, is_climate := climate_relevance ==1L]godad[, is_nonclimate := climate_relevance ==0L]godad[, is_mitigation :=str_detect(meta_category, "\\bmitig")] # robustgodad[, is_adaptation :=str_detect(meta_category, "\\badapt")] # robust# helper: first adoption year by ADM2first_adopt <-function(dt_subset) {if (nrow(dt_subset) ==0L) return(data.table(adm2_id =character(), g =integer())) out <- dt_subset[, .(g =suppressWarnings(min(as.integer(year), na.rm =TRUE))), by = adm2_id] out[is.infinite(g), g :=NA_integer_][]}adopt_adapt <-first_adopt(godad[is_adaptation ==TRUE])adopt_mitig <-first_adopt(godad[is_mitigation ==TRUE])adopt_climate <-first_adopt(godad[is_climate ==TRUE])adopt_nonclimate <-first_adopt(godad[is_nonclimate ==TRUE])sapply(list(adapt = adopt_adapt, mitig = adopt_mitig, climate = adopt_climate, nonclimate = adopt_nonclimate), nrow)``````{r}#| label: csdid-wrapperyvar <-if ("pollution_CO2_log"%in%names(base_panel)) "pollution_CO2_log"else"pollution_CO2_w"run_csdid <-function(base_panel, adopt_tbl, label, min_e =-10, max_e =20) { P <-merge(base_panel[, .(adm2_id, year, y =get(yvar))], adopt_tbl, by ="adm2_id", all.x =TRUE)# Treated flags & event time P[, treated :=as.integer(!is.na(g) & year >= g)] P[, rel_time :=ifelse(is.na(g), NA_integer_, year - g)] P[, adm2_id_int :=as.integer(factor(adm2_id))]# Drop missing outcome and units treated in first overall year min_year <-min(P$year, na.rm =TRUE) D <- P[!is.na(y) & (is.na(g) | g > min_year)]if (nrow(D[!is.na(g)]) ==0L) {warning(sprintf("No treated units for %s — skipping.", label))return(NULL) } att <- did::att_gt(yname ="y",tname ="year",idname ="adm2_id_int",gname ="g",data = D,panel =TRUE, # repeated cross-sections at ADM2-yearcontrol_group ="notyettreated",allow_unbalanced_panel =TRUE )list(label = label,att = att,es = did::aggte(att, type ="dynamic", min_e = min_e, max_e = max_e),grp = did::aggte(att, type ="group") )}```# Empirical Strategy## Staggered adoption DiD (Callaway–Sant’Anna)Let $Y_{jt}$ denote log pollution in ADM2 $j$ and year $t$. We estimate group-time average treatment effects using the **CS** estimator with not-yet-treated controls:$ATT(g,e) = \mathbb{E}[Y_{jt}(1) - Y_{jt}(0) \mid g_j = g,, t-g = e],$and aggregate to dynamic event-time profiles with simultaneous confidence bands.### Treatment intensityTo study **heterogeneity by dose**, we compute **post-treatment totals** (commitments/disbursements) and re-estimate dynamics **within dose bins**. This recovers a policy-relevant “more money means larger effects?” gradient while preserving the staggered-adoption identification. As a complementary specification, we model **continuous doses** in a two-way fixed-effects panel with distributed lags of the annual amount in $j,t$.# Sample description ```{r}#| label: descriptives-setup# Build a descriptive-sample panel using "All climate" adoption (change if desired)descP <-merge( base_panel[, .(adm2_id, year, pollution_CO2, pollution_CO2_log)], adopt_climate, by ="adm2_id", all.x =TRUE)descP[, treated :=as.integer(!is.na(g) & year >= g)]descP[, ever_treated :=as.integer(!is.na(g))]# Helper: nice quantile summaryq_summ <-function(x) {c(N =sum(!is.na(x)),mean =mean(x, na.rm =TRUE),sd =sd(x, na.rm =TRUE),p10 =quantile(x, .10, na.rm =TRUE),p50 =quantile(x, .50, na.rm =TRUE),p90 =quantile(x, .90, na.rm =TRUE))}```## Sample overview (counts, years, treated share)```{r}#| label: descriptives-overviewlibrary(data.table)ov <-list(N_obs =nrow(descP),N_ids =uniqueN(descP$adm2_id),N_years =uniqueN(descP$year),year_min_max =paste(range(descP$year, na.rm =TRUE), collapse ="–"),N_ever_treated = descP[, uniqueN(adm2_id[ever_treated ==1])],N_never_treated = descP[, uniqueN(adm2_id[ever_treated ==0|is.na(ever_treated)])],share_treated_by_2023 =with(descP[year ==max(year, na.rm =TRUE)],mean(ever_treated ==1, na.rm =TRUE)))as.data.table(ov)```## Outcome summaries (overall & by ever-treated)```{r}#| label: descriptives-outcome-summariessumm_all <-q_summ(descP$pollution_CO2_log)by_status <- descP[, as.list(q_summ(pollution_CO2_log)), by = ever_treated]as.data.table(t(summ_all))[]by_status[]```## Adoption timing (cohort sizes, cumulative share)```{r}#| label: descriptives-adoptionlibrary(ggplot2)# Cohort histogram (first g per ADM2)cohort_sizes <- adopt_climate[!is.na(g), .N, by = g][order(g)]cohort_sizes%>%filter(g>=2000)%>%ggplot(aes(g, N)) +geom_col() +labs(title ="Cohort size by first treatment year (All climate)",x ="First project year (g)", y ="ADM2 count") +theme_classic(base_size =12)# Cumulative treated share over time# Cumulative treated share = fraction of ADM2 whose g <= yearshare_treated <- descP[!is.na(g), .(share =mean(g <= year)), by = year]# Also include never-treated in denominatorall_ids <-unique(descP$adm2_id)total_n <-length(all_ids)share_treated <- descP[, .(share =mean(!is.na(g) & g <= year)), by = year]ggplot(share_treated, aes(year, share)) +geom_line() +geom_point() +scale_y_continuous(labels = scales::percent) +labs(title ="Share of ADM2 ever treated (cumulative)",x =NULL, y =NULL) +theme_classic(base_size =12)```## Outcome distributions & time trends```{r}# Density of log outcome overallggplot(descP[!is.na(pollution_CO2_log)], aes(pollution_CO2_log)) +geom_density() +labs(title ="Distribution of log(1+CO₂) across ADM2-years",x ="pollution_CO2_log", y ="Density") +theme_classic(base_size =12)# Density by ever-treated statusggplot(descP[!is.na(pollution_CO2_log)],aes(pollution_CO2_log, fill =factor(ever_treated))) +geom_density(alpha =0.35) +scale_fill_discrete(name ="Ever treated") +labs(title ="Distribution of log(1+CO₂) by ever-treated status",x ="pollution_CO2_log", y ="Density") +theme_classic(base_size =12)# Mean log outcome over time by ever-treated statusmean_by_year <- descP[!is.na(pollution_CO2_log), .(y =mean(pollution_CO2_log, na.rm =TRUE)), by = .(year, ever_treated)]ggplot(mean_by_year, aes(year, y, color =factor(ever_treated))) +geom_line() +geom_point(size =0.9) +labs(title ="Mean log(1+CO₂) over time",color ="Ever treated", x =NULL, y ="Mean log(1+CO₂)") +theme_classic(base_size =12)```## Raw event-time mean (not causal; descriptive)```{r}#| label: descriptives-event-time-rawdescP[, rel_time :=ifelse(is.na(g), NA_integer_, year - g)]raw_es <- descP[!is.na(rel_time), .(y =mean(pollution_CO2_log, na.rm =TRUE)), by = rel_time]ggplot(raw_es, aes(rel_time, y)) +geom_line() +geom_point() +geom_vline(xintercept =-1, linetype =2) +labs(title ="Raw event-time mean: log(1+CO₂) (All climate)",x ="Event time (t - g)", y ="Mean log(1+CO₂)") +theme_classic(base_size =12)```{fig-align="center"}# Main dynamic effects (overall event study)We identify the **first year** in which an ADM2 receives at least one **climate** project. If the GODAD file already has an ADM2 code, we use it; if not, we perform a spatial join using the point geometry.```{r}#| label: csdid-runres_adapt <-run_csdid(base_panel, adopt_adapt, "Adaptation")res_mitig <-run_csdid(base_panel, adopt_mitig, "Mitigation")res_clim <-run_csdid(base_panel, adopt_climate, "All climate")res_noncl <-run_csdid(base_panel, adopt_nonclimate, "Non-climate")``````{r}#| label: csdid-plot-sum#| fig.width: 12#| fig.height: 10make_es_df <-function(res) {if (is.null(res)) return(NULL) es <- res$esdata.frame(label = res$label, egt = es$egt, att = es$att.egt, se = es$se.egt)}es_all <-rbind(make_es_df(res_adapt),make_es_df(res_mitig),make_es_df(res_clim),make_es_df(res_noncl))if (!is.null(es_all)) {ggplot(es_all, aes(egt, att)) +geom_hline(yintercept =0, linetype =2) +geom_point() +geom_errorbar(aes(ymin = att -1.96*se, ymax = att +1.96*se), width =0.15) +geom_vline(xintercept =-1, linetype =2) +facet_wrap(~ label, ncol =2, scales ="free_y") +labs(title =sprintf("Dynamic ATT (y = %s)", yvar),x ="Event time (t - g)", y ="ATT") +theme_minimal(base_size =12)}summ_line <-function(res) {if (is.null(res)) return(NULL) s <-capture.output(print(summary(res$grp)))data.frame(model = res$label, summary =paste(s, collapse ="\n"))}summ_all <- dplyr::bind_rows(summ_line(res_adapt),summ_line(res_mitig),summ_line(res_clim),summ_line(res_noncl))if (nrow(summ_all)) {cat("\n\n===== Overall ATT (group aggregated) =====\n")for (i inseq_len(nrow(summ_all))) {cat("\n--", summ_all$model[i], "--\n", summ_all$summary[i], "\n") }}# SAFE overall row from a result objectsumm_line_overall <-function(res, use_fallback =TRUE) {if (is.null(res) ||is.null(res$grp)) return(NULL)# Prefer group overall directly without calling summary() (avoids noisy prints) att <-tryCatch(res$grp$overall.att, error =function(e) NA_real_) se <-tryCatch(res$grp$overall.se, error =function(e) NA_real_)if (length(att) ==1L &&is.finite(att) &&length(se) ==1L &&is.finite(se)) {return(data.frame(label = res$label, att = att, se = se)) }# Optional fallback: inverse-variance weighted mean of post dynamic ATTsif (use_fallback &&!is.null(res$es) &&length(res$es$att.egt)) { df <-data.frame(att = res$es$att.egt, se = res$es$se.egt, e = res$es$egt) df <- df[is.finite(df$att) &is.finite(df$se) & df$e >=0, ]if (nrow(df) >=1&&all(df$se >0)) { w <-1/ (df$se^2) att <-sum(w * df$att) /sum(w) se <-sqrt(1/sum(w))return(data.frame(label = res$label, att = att, se = se)) } }# Nothing usable for this specNULL}summ_all <- dplyr::bind_rows(summ_line_overall(res_adapt), summ_line_overall(res_mitig),summ_line_overall(res_clim), summ_line_overall(res_noncl)) |> dplyr::mutate(pct =bt(att),lo =bt(att, se)["lo"], hi =bt(att, se)["hi"] )gt::gt(summ_all) |> gt::fmt_number(columns =c(pct, lo, hi), decimals =1) |> gt::cols_label(pct="ATT (%)", lo="95% lo", hi="95% hi") |> gt::tab_caption("Overall ATT (group aggregated), back-transformed to %.")```**Notes.** The absence of positive pre-trends (e \< 0) supports the parallel-trends assumption.# Intensity: event study by post-treatment dose bins```{r}#| label: amounts-from-godad-setup#| message: false#| warning: falsesuppressPackageStartupMessages({library(data.table); library(ggplot2); library(did)})setDT(godad); setDT(base_panel)# Prefer per-location split columns if already computed in `godad`amt_cols_priority <-c("disb_loc_evensplit", "comm_loc_evensplit","disb_amount", "commit_amount", "amount")amount_col_godad <- amt_cols_priority[amt_cols_priority %in%names(godad)][1]stopifnot(length(amount_col_godad) ==1)message("Using amount column from godad: ", amount_col_godad)# Sanity: we need ADM2 id and YEAR in godad; if year is a date, coerce to integerif (!"year"%in%names(godad)) {stop("`godad` must have a `year` column (integer).")}if (!"adm2_id"%in%names(godad)) {stop("`godad` must carry `adm2_id` (or do the spatial join before this step).")}# Helper to build ADM2-year totals for a filteradm2yr_sum <-function(gd, keep_expr) { x <- gd[eval(keep_expr), .(amt =sum(get(amount_col_godad), na.rm =TRUE)), by = .(adm2_id, year)]setnames(x, "amt", "dose_amt") x[]}# Category filters (adjust if your flags differ)gd <-copy(godad)has_cols <-names(gd)req_flags <-c("is_adaptation","is_mitigation","is_climate")if (!all(req_flags %in% has_cols)) {stop("godad must have logical flags: is_adaptation, is_mitigation, is_climate.")}adm2yr_adapt <-adm2yr_sum(gd, quote(is_adaptation ==TRUE))adm2yr_mitig <-adm2yr_sum(gd, quote(is_mitigation ==TRUE))adm2yr_climate <-adm2yr_sum(gd, quote(is_climate ==TRUE))adm2yr_nonclimate <-adm2yr_sum(gd, quote(is_climate ==FALSE))# Generic merge + first adoption helperattach_dose_and_adopt <-function(bp, adm2yr) { dt <- adm2yr[bp, on =c("adm2_id","year")] dt[is.na(dose_amt), dose_amt :=0] # zero when no project that year# First year with positive amount is the adoption g gtab <- dt[dose_amt >0, .(g =min(year)), by = adm2_id] dt <- gtab[dt, on ="adm2_id"] dt[]}bp_adapt <-attach_dose_and_adopt(base_panel, adm2yr_adapt)bp_mitig <-attach_dose_and_adopt(base_panel, adm2yr_mitig)bp_climate <-attach_dose_and_adopt(base_panel, adm2yr_climate)bp_nonclimate <-attach_dose_and_adopt(base_panel, adm2yr_nonclimate)post_totals_and_bins <-function(dt, n_bins =4) { d <-copy(dt)setorder(d, adm2_id, year)# cumulative dose (for convenience) d[, dose_cum :=cumsum(dose_amt), by = adm2_id]# total post-treatment dose per unit (from g onward) posttab <- d[!is.na(g), .(post_total =max(dose_cum[year >= g[1L]], na.rm =TRUE)), by = adm2_id] d <- posttab[d, on ="adm2_id"]# Cut bins on treated units only treated <- d[!is.na(g), unique(adm2_id)] x <- d[adm2_id %in% treated, post_total]# robust breaks mk_breaks <-function(x, n_bins =4) {if (length(unique(na.omit(x))) < n_bins) x <- x +rnorm(length(x), sd =sd(x, na.rm=TRUE)*1e-8) b <-quantile(x, probs =seq(0,1,length.out = n_bins+1), na.rm =TRUE, type =1) |>unique() |>sort()if (length(b) <=2) b <-quantile(x, probs =c(0,.5,1), na.rm =TRUE, type =1) |>unique() |>sort() b } brks <-mk_breaks(x, n_bins) d[, post_bin :=cut(post_total, breaks = brks, include.lowest =TRUE, right =TRUE)] d[, post_bin :=droplevels(post_bin)] d[]}bp_adapt <-post_totals_and_bins(bp_adapt, n_bins =4)bp_mitig <-post_totals_and_bins(bp_mitig, n_bins =4)bp_climate <-post_totals_and_bins(bp_climate, n_bins =4)bp_nonclimate <-post_totals_and_bins(bp_nonclimate, n_bins =4)``````{r}#| label: csdid-postdose-module#| message: false#| warning: false#| fig-cap: "Event-study by post-treatment total dose — stratified bins (treated-only cuts)."yvar <-"pollution_CO2_log"; stopifnot(yvar %in%names(base_panel))suppressPackageStartupMessages({ library(data.table); library(did) })check_and_fix_postdose <-function(dtf, yvar) { d <-as.data.table(copy(dtf))# 1) Outcome present?if (!yvar %in%names(d)) stop("Outcome column `", yvar, "` not found in data.")# 2) Year must be integerif (!is.integer(d$year)) { d[, year :=as.integer(year)]if (anyNA(d$year)) stop("`year` coercion produced NAs — check your year values.") }# 3) ID should be numeric for did stabilityif (!is.numeric(d$adm2_id)) { d[, adm2_id_int :=as.integer(factor(adm2_id))] } else {setnames(d, "adm2_id", "adm2_id_int") }# 4) Event-time origin g (first treat year) must exist for treatedif (!"g"%in%names(d)) stop("`g` (first treatment year) is missing.")if (all(is.na(d$g))) warning("All g are NA — this category may have no treated units.")# 5) post_bin must exist and have levels with treated obsif (!"post_bin"%in%names(d)) stop("`post_bin` is missing. Run the bin construction step first.") d[, post_bin :=droplevels(post_bin)]if (nlevels(d$post_bin) ==0L) stop("All `post_bin` are NA — no treated ADM2s or post_total is missing.")# Keep only bins that actually contain treated units has_tr <- d[, any(!is.na(g)), by = post_bin] valid_bins <- has_tr[ V1 ==TRUE, post_bin ] d <- d[post_bin %in% valid_bins] d[, post_bin :=droplevels(post_bin)]# 6) Drop rows with missing outcome d <- d[!is.na(get(yvar))]# 7) Basic summary to the consoleprint(d[, .(n_rows=.N,n_adm2=uniqueN(adm2_id_int),n_treated=uniqueN(adm2_id_int[!is.na(g)]),bins=nlevels(post_bin))])print(d[, .(n_adm2=uniqueN(adm2_id_int),n_treated=uniqueN(adm2_id_int[!is.na(g)])),by=post_bin][order(post_bin)])return(d[])}# Wrapper to run CS-DiD with the numeric id we just ensuredrun_postdose_csdid_safe <-function(dtf, label, yvar, min_e=-10, max_e=20) { dd <-check_and_fix_postdose(dtf, yvar) out_es <-vector("list", length(levels(dd$post_bin)))names(out_es) <-levels(dd$post_bin)for (b inlevels(dd$post_bin)) { sub <- dd[post_bin == b]if (sub[!is.na(g), .N] ==0L) next att <- did::att_gt(yname = yvar,tname ="year",idname ="adm2_id_int", # <- numeric idgname ="g",data = sub,panel =TRUE,control_group ="notyettreated",bstrap =TRUE,clustervars ="adm2_id_int" ) es <- did::aggte(att, type ="dynamic", min_e = min_e, max_e = max_e, na.rm = T) out_es[[b]] <- es } es_df <- data.table::rbindlist(lapply(names(out_es), function(b){ x <- out_es[[b]]; if (is.null(x)) return(NULL)data.table(bin=b, e=x$egt, att=x$att.egt, se=x$se.egt) }), use.names=TRUE, fill=TRUE) bin_order <- dd[, .(ord =median(post_total, na.rm =TRUE)), by = post_bin][order(ord), as.character(post_bin)] es_df[, bin :=factor(bin, levels = bin_order, ordered =TRUE)] # <<< dd[, post_bin :=factor(post_bin, levels = bin_order, ordered =TRUE)] # <<< (keeps summaries consistent)if (nrow(es_df)) {library(ggplot2)print(ggplot(es_df, aes(x = e, y = att)) +# horizontal line at 0geom_hline(yintercept =0, linetype =2) +# error barsgeom_errorbar(aes(ymin = att -1.96* se, ymax = att +1.96* se),width =0.15) +# pointsgeom_point() +# vertical line at -1 (pre-treatment period marker)geom_vline(xintercept =-1, linetype =2) +# facets by intensity binfacet_wrap(~ bin, scales ="free_y") +labs(x ="Event time (years since first project)",y ="ATT on outcome",title =paste("CS-DiD dynamics by post-treatment dose —", label) ) +theme_minimal(base_size =12) +theme(strip.text =element_text(face ="bold"),plot.title =element_text(face ="bold"),panel.grid.minor =element_blank() ) ) }invisible(list(es=out_es, df=es_df))}# Run all four familiesres_adapt <-run_postdose_csdid_safe(bp_adapt, "Adaptation", yvar)res_mitig <-run_postdose_csdid_safe(bp_mitig, "Mitigation", yvar)res_climate <-run_postdose_csdid_safe(bp_climate, "All climate", yvar)res_nonclimate <-run_postdose_csdid_safe(bp_nonclimate, "Non-climate", yvar)```**Interpretation.** Higher post-treatment exposure is associated with larger post-event estimates, consistent with a **dose–response** relationship. Pre-treatment coefficients remain centered at zero across bins.## Excluding China```{r}#| message: false#| warning: falseadopt_adapt <- adopt_adapt[!grepl("CHN", adm2_id), ]adopt_mitig <- adopt_mitig[!grepl("CHN", adm2_id), ]adopt_climate <- adopt_climate[!grepl("CHN", adm2_id), ]adopt_nonclimate <- adopt_nonclimate[!grepl("CHN", adm2_id), ]res_adapt_china <-run_csdid(base_panel, adopt_adapt, "Adaptation")res_mitig_china <-run_csdid(base_panel, adopt_mitig, "Mitigation")res_clim_china <-run_csdid(base_panel, adopt_climate, "All climate")res_noncl_china <-run_csdid(base_panel, adopt_nonclimate, "Non-climate")es_all_china <-rbind(make_es_df(res_adapt_china),make_es_df(res_mitig_china),make_es_df(res_clim_china),make_es_df(res_noncl_china))if (!is.null(es_all_china)) {ggplot(es_all_china, aes(egt, att)) +geom_hline(yintercept =0, linetype =2) +geom_point() +geom_errorbar(aes(ymin = att -1.96*se, ymax = att +1.96*se), width =0.15) +geom_vline(xintercept =-1, linetype =2) +facet_wrap(~ label, ncol =2, scales ="free_y") +labs(title =sprintf("Dynamic ATT (y = %s)", yvar),x ="Event time (t - g)", y ="ATT") +theme_minimal(base_size =12)}summ_all_china <- dplyr::bind_rows(summ_line(res_adapt_china),summ_line(res_mitig_china),summ_line(res_clim_china),summ_line(res_noncl_china))if (nrow(summ_all)) {cat("\n\n===== Overall ATT (group aggregated) =====\n")for (i inseq_len(nrow(summ_all_china))) {cat("\n--", summ_all_china$model[i], "--\n", summ_all_china$summary[i], "\n") }}summ_all_china <- dplyr::bind_rows(summ_line_overall(res_adapt), summ_line_overall(res_mitig),summ_line_overall(res_clim), summ_line_overall(res_noncl)) |> dplyr::mutate(pct =bt(att),lo =bt(att, se)["lo"], hi =bt(att, se)["hi"] )gt::gt(summ_all_china) |> gt::fmt_number(columns =c(pct, lo, hi), decimals =1) |> gt::cols_label(pct="ATT (%)", lo="95% lo", hi="95% hi") |> gt::tab_caption("Overall ATT (group aggregated), back-transformed to %.")```## Intensity Adapt excluding China```{r}#| label: intensity-dose-response#| echo: false#| message: false#| warning: false# Run all four familiesres_adapt <-run_postdose_csdid_safe(bp_adapt[!grepl("CHN", adm2_id), ], "Adaptation", yvar)res_mitig <-run_postdose_csdid_safe(bp_mitig[!grepl("CHN", adm2_id), ], "Mitigation", yvar)res_climate <-run_postdose_csdid_safe(bp_climate[!grepl("CHN", adm2_id), ], "All climate", yvar)res_nonclimate <-run_postdose_csdid_safe(bp_nonclimate[!grepl("CHN", adm2_id), ], "Non-climate", yvar)``````{r}#| label: descriptives-dists-trends# Density of log outcome overallggplot(descP[!is.na(pollution_CO2_log)], aes(pollution_CO2_log)) +geom_density() +labs(title ="Distribution of log(1+CO₂) across ADM2-years",x ="pollution_CO2_log", y ="Density") +theme_classic(base_size =12)# Density by ever-treated statusggplot(descP[!is.na(pollution_CO2_log)],aes(pollution_CO2_log, fill =factor(ever_treated))) +geom_density(alpha =0.35) +scale_fill_discrete(name ="Ever treated") +labs(title ="Distribution of log(1+CO₂) by ever-treated status",x ="pollution_CO2_log", y ="Density") +theme_classic(base_size =12)# Mean log outcome over time by ever-treated statusmean_by_year <- descP[!is.na(pollution_CO2_log), .(y =mean(pollution_CO2_log, na.rm =TRUE)), by = .(year, ever_treated)]ggplot(mean_by_year, aes(year, y, color =factor(ever_treated))) +geom_line() +geom_point(size =0.9) +labs(title ="Mean log(1+CO₂) over time",color ="Ever treated", x =NULL, y ="Mean log(1+CO₂)") +theme_classic(base_size =12)``````{# {r}# #| label: map-posttotal# #| fig-cap: "Total post-treatment adaptation amounts by ADM2"# #| warning: false# #| message: false# # # library(sf)# library(ggplot2)# library(scales)# # adm2$adm2_id <- adm2$GID_2# # map_df <- adm2 |># merge(bp_adapt[, .(post_total = max(post_total, na.rm = TRUE)), by = adm2_id],# by = "adm2_id", all.x = TRUE)# # world_bg <- rnaturalearth::ne_countries(scale = "medium", returnclass = "sf") |># st_transform(st_crs(map_df)) # make sure CRS matches your ADM2 layer# # map <- ggplot() +# # Grey background countries# geom_sf(data = world_bg, fill = "grey90", color = "white", size = 0.1) +# # # Your ADM2 layer# geom_sf(data = map_df, aes(fill = post_total), color = NA) + # scale_fill_viridis_c(# option = "C",# trans = "log",# labels = label_number(scale_cut = cut_si("unit")),# na.value = "grey90",# name = "Total\n(Climate Finance)"# ) +# labs(# title = "Where is post-treatment Climate finance concentrated?",# subtitle = "ADM2-level total amounts, log scale. Light grey = no projects.",# caption = "Source: GODAD; Author's calculations."# ) +# theme_minimal(base_size = 11) +# theme(# panel.grid.major = element_blank(),# panel.grid.minor = element_blank(),# legend.position = "right",# plot.title = element_text(face = "bold"),# plot.subtitle = element_text(size = 9)# )# # ggsave("Figures/Climate_adm2.png", plot = map, width = 20, height = 15.22, units = "in", dpi = 1000)``````{# {r}# #| label: waffle-dose-bins# #| fig-cap: "Share of total adaptation amount by dose bin (1 square = ~1M)"# #| warning: false# #| message: false# library(waffle)# library(dplyr)# library(stringr)# library(RColorBrewer)# # # 1️⃣ Aggregate total amounts by dose bin# bin_amounts <- bp_adapt[!is.na(post_bin),# .(total_amt = sum(dose_amt, na.rm = TRUE)), # by = post_bin] %>%# mutate(post_bin = as.character(post_bin))# # # 2️⃣ Extract numeric lower bound from bin labels for ordering# get_lower <- function(x) {# as.numeric(str_extract(x, "[-+]?[0-9]*\\.?[0-9]+(?:e[+-]?\\d+)?"))# }# # bin_amounts <- bin_amounts %>%# mutate(lower_bound = get_lower(post_bin)) %>%# arrange(lower_bound)# # # 3️⃣ Convert to millions and build named vector in the **sorted order**# parts <- round(bin_amounts$total_amt / 1e7) # 1 square = 10M# names(parts) <- bin_amounts$post_bin# # # 4️⃣ Define number of bins explicitly# n_bins <- length(parts)# # # 5️⃣ Pick a nice palette with the right length# # Paired works well up to 12; or use viridis(n_bins)# pal <- brewer.pal(n_bins, "Paired") # or "PuBuGn", "YlOrRd", "Set2"...# # bin_levels <- names(parts) # order from lowest to highest bin# # waffle_df <- tibble(# part = factor(rep(bin_levels, times = parts), levels = bin_levels)# ) %>%# mutate(# index = row_number(),# y = (index - 1) %/% 20, # number of rows# x = (index - 1) %% 20# )# # # Reverse y so origin is bottom-left# waffle_df$y <- max(waffle_df$y) - waffle_df$y# # # --- Plot ---# waffle <- ggplot(waffle_df, aes(x, y, fill = part)) +# geom_tile(color = "white", size = 0.25) +# coord_equal() +# scale_fill_manual(values = pal, name = "Dose bin") +# labs(# title = "Share of total adaptation amount by dose bin",# subtitle = "Each square represents approximately 10 million USD.\nColors represent bins of post-treatment finance intensity (lower → higher).",# x = NULL, y = NULL,# caption = "Source: GODAD; Author's calculations."# ) +# theme_minimal(base_size = 12) +# theme(# axis.text = element_blank(),# axis.ticks = element_blank(),# panel.grid = element_blank(),# legend.position = "right",# plot.title = element_text(face = "bold"),# plot.subtitle = element_text(size = 10),# plot.caption = element_text(size = 8, hjust = 0)# )# # ggsave(# "Figures/Adaptation_waffle.png",# dpi=1000, width = 12, height = 9# )``````{r}library(data.table)setDT(godad)# --- Aggregate by country-year ---country_year_totals <- godad[, .(adapt_total =sum(comm_loc_evensplit[is_adaptation], na.rm =TRUE),mitig_total =sum(comm_loc_evensplit[is_mitigation], na.rm =TRUE),climate_total =sum(comm_loc_evensplit[is_climate], na.rm =TRUE),nonclimate_total =sum(comm_loc_evensplit[is_nonclimate], na.rm =TRUE)), by = .(gid_0, year)]```\newpage# References {.unnumbered}::: {#refs}:::\newpage# Appendix {.unnumbered}\setcounter{table}{0}\renewcommand{\thetable}{A\arabic{table}}\setcounter{figure}{0}\renewcommand{\thefigure}{A\arabic{figure}}## Top emitters (ADM2) — average over the sample window```{r}#| label: descriptives-top-emitters# Extract ISO3 from GADM code pattern like "CHN.10.9_1"descP[, iso3 :=sub("\\..*$", "", adm2_id)]top_adm2 <- descP[, .(avg_CO2 =mean(pollution_CO2, na.rm =TRUE)), by = .(iso3, adm2_id)][order(-avg_CO2)][1:20]top_adm2[]```## Heatmaps (country × year matrix)```{r}library(forcats)country_year_heat <- country_year_totals %>%mutate(country =fct_reorder(gid_0, adapt_total, .fun = max, .desc =TRUE))ggplot(country_year_heat, aes(x = year, y = country, fill =log1p(adapt_total))) +geom_tile() +scale_fill_viridis_c() +labs(fill ="Log Adaptation\nFinance")```{fig-align="center"}## Reproducibility```{r}#| label: session-info#| echo: falsesessionInfo()``````{r}library(ggridges)ggplot(country_year_totals, aes(x =log1p(nonclimate_total), y =factor(year))) +geom_density_ridges(scale =3, rel_min_height =0.01) +labs(x ="Log Adaptation Finance", y ="Year")``````{r}library(dplyr)library(tidyr)library(ggplot2)library(ggbump)# rank each yearranked <- country_year_totals %>%filter(year %in%1995:2023)%>%mutate(year =as.integer(year)) %>%group_by(year) %>%mutate(rank_y =min_rank(desc(adapt_total))) %>%ungroup()# contenders = countries that appear in top10 at least 3 yearscontenders <- ranked %>%filter(rank_y <=10) %>%count(gid_0, name ="n_top10") %>%filter(n_top10 >=3) %>%pull(gid_0)# build continuous panel for contenders (fill missing years with 0)years_all <-sort(unique(ranked$year))panel <- country_year_totals %>%filter(year %in%1995:2023)%>%filter(gid_0 %in% contenders) %>%select(gid_0, year, adapt_total) %>%complete(gid_0, year = years_all, fill =list(adapt_total =0)) %>%group_by(year) %>%mutate(rank_c =min_rank(desc(adapt_total))) %>%# rank among contendersungroup()ggplot(panel, aes(year, rank_c, color = gid_0, group = gid_0)) +geom_bump(size =1.2) +geom_point(size =1.8) +scale_y_reverse(breaks =1:10) +labs(title ="Top adaptation finance recipients over time (contenders)",x =NULL, y ="Rank") +theme_minimal(base_size =12) +theme(legend.position ="none",plot.title =element_text(face ="bold"))``````{r}library(ggbump)# build ranks per yearranked <- country_year_totals %>%mutate(year =as.integer(year)) %>%group_by(year) %>%mutate(rank_y =min_rank(desc(adapt_total))) %>%ungroup()# take top 10 each yeartop10 <- ranked %>%filter(rank_y <=10)# those with at least 2 years in top10 (for smooth bumps)lines_df <- top10 %>%group_by(gid_0) %>%filter(n() >=2) %>%ungroup()ggplot(lines_df, aes(x = year, y = rank_y, color = gid_0, group = gid_0)) +geom_bump(size =1.2) +geom_point(size =2) +scale_y_reverse(breaks =1:10) +labs(title ="Top 10 adaptation finance recipients over time",x =NULL, y ="Rank",color ="Country" ) +theme_minimal(base_size =12) +theme(plot.title =element_text(face ="bold"),legend.position ="right" )```
